
### Doubly robust estimating equation
## Input: 
##      current estimate of beta (params for logop model) and prop score
##      va:     covariates for log rr model; vb: covariates for log op model
##      pscore: a vector giving propensity scores for each observation
##      wt:     is a set of weights that can be set for optimal efficiency
## Output: 
##      New estimate of alpha (params for logrr model)


dr.objective.iter <- function(par.update, parn = "not", pars, data, cnt){
        # pars represents alpha here. beta is beta, gamma is gamma.
        # startpars <- 
        
        # pars <- startpars <- alpha.true
        # beta <- beta.true; gamma <- gamma.true
        l0 = data[,c("l0.1","l0.2")]
        a0 = data[,"a0"];  l1 = data[,"l1"];  a1 = data[,"a1"];   y = data[,"y"]
        nb = nrow(data)
        
        if(parn == "alpha1")     pars[1,]   = par.update
        if(parn == "alpha2345")  pars[2:5,] = par.update
        
        base.covs <- func.DataGen(c(pars, c(beta), c(gamma)), data,cnt)
        MyAttach(base.covs)
        
        theta.k <- sapply(1:nb, function(i) theta[i,k[i]])
        theta.k.inv <- 1/theta.k
        theta.k.inv.a1 <- (theta.k.inv)^a1
        # for (i in 1:nb){theta.k[i] <- theta[i,][k[i]]}
        
        ### m = 1
        u1 <- y * theta.k.inv.a1 
        E.u1 <- g(l1, a0, l0, p.y.cond, p.a1, theta.k.inv, gamma, nb)
        
        d1.row.k <- -l0 * a1 * p.y * theta.k.inv.a1
        d1 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
        for (i in 1:nb) {
            d1[i,k[i],] <- d1.row.k[i,]
        }
        
        
        # E.d1.row.k <- -l0 * p.a1 * E.u1
        E.d1.row.k <- -l0 * g1(l1, a0, l0, p.y.cond, p.a1, theta.k.inv, gamma, nb)
        E.d1 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
        for (i in 1:nb) {
            E.d1[i,k[i],] <- E.d1.row.k[i,]
        }
        # E.d1[1,,]
        
        
        ### m = 0
        theta. <- sapply(1:nb, function(i) theta[i,1])
        theta..inv <- 1/theta.
        theta..inv.a0 <- (theta..inv)^a0
        u0 <- y * theta.k.inv.a1 * theta..inv.a0
        E.u0 <- p.a0 * h(a0=vec.1, l0, theta..inv.a0, 
                         p.y.cond, p.a1, theta.k.inv, gamma, nb,
                         phi45) + 
            (1-p.a0) * h(a0=vec.0, l0, theta..inv.a0, 
                         p.y.cond, p.a1, theta.k.inv, gamma, nb,
                         phi45)
        
        d0.row.1 <- f.(a0,l0, theta..inv.a0, 
                       p.y.cond, p.a1, theta.k.inv, gamma, nb,
                       phi45)
        d0.row.k <- fk(a0,l0, theta..inv.a0, 
                       p.y.cond, p.a1, theta.k.inv, gamma, nb,
                       phi45)
        d0 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
        for (i in 1:nb) {
            d0[i,1,] <- d0.row.1[i,]
            d0[i,k[i],] <- d0.row.k[i,]
        }
        
        E.d0.row.1 <- p.a0 * f.(a0=vec.1,l0, theta..inv.a0, 
                                p.y.cond, p.a1, theta.k.inv, gamma, nb,
                                phi45) +
            (1- p.a0) * f.(a0=vec.0,l0, theta..inv.a0, 
                           p.y.cond, p.a1, theta.k.inv, gamma, nb,
                           phi45)
        E.d0.row.k <- p.a0 * fk(a0=vec.1,l0, theta..inv.a0, 
                                p.y.cond, p.a1, theta.k.inv, gamma, nb,
                                phi45) +
            (1- p.a0) * fk(a0=vec.0,l0, theta..inv.a0, 
                           p.y.cond, p.a1, theta.k.inv, gamma, nb,
                           phi45)
        E.d0 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
        for (i in 1:nb) {
            E.d0[i,1,] <- E.d0.row.1[i,]
            E.d0[i,k[i],] <- E.d0.row.k[i,]
        }
        
        
        tmp <-  array(rep(0,5*2),dim=c(5, 2))
        tmp1 <-  array(rep(0,5*2),dim=c(5, 2)) # DEBUG
        tmp3 <- 0
        # tmp2 <- matrix(rep(0,32*2),ncol=2) # 
        for (i in 1:nb){
            # total <- (d1 - E.d1) * (u1 - E.u1) + (d0 - E.d0) * (u0 - E.u0) 
            total <- (d1 - E.d1) * (u1 - E.u1)
            tmp <- tmp + (total[i,,])^2 * cnt[i]
            tmp1 <- tmp1 + total[i,,] * cnt[i]
            tmp3 <- tmp3 + (u0 - E.u0)[i]
            # tmp2 <- tmp2 + (d0.row.1  - E.d0.row.1 )*cnt[i]
        }
        return(sum(tmp))
}
